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Abstract 

The speed of many one-line transformation methods for the produc- 
tion of, for example, Levy alpha-stable random numbers, which generalize 
Gaussian ones, and Mittag-LefHer random numbers, which generalize ex- 
ponential ones, is very high and satisfactory for most purposes. However, 
for the class of decreasing probability densities fast rejection implemen- 
tations like the Ziggurat by Marsaglia and Tsang promise a significant 
speed-up if it is possible to complement them with a method that sam- 
ples the tails of the infinite support. This requires the fast generation of 
random numbers greater or smaller than a certain value. We present a 
method to achieve this, and also to generate random numbers within any 
arbitrary interval. We demonstrate the method showing the properties 
of the transform maps of the above mentioned distributions as examples 
of stable and geometric stable random numbers used for the stochastic 
solution of the space-time fractional diffusion equation. 

February 18, 2009 



1 Introduction 

Many numerical methods for the generation of random numbers represent the 
main body of the probability density using a fast method and the tails using an 
alternative method. A famous example is the Ziggurat method by Marsaglia and 
Tsang [57]. Fig. fl] depicts the situation schematically. A reason for this apparent 
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Figure 1 : Schematic illustration of using two methods for sampling a distribu- 
tion. A fast method may be available only for the body while the tails can be 
sampled with a slower method. 

complication is that the method for the main body works best and fastest on a 
finite support or is specially designed for the main body in terms of accuracy or 
speed. Handling the tails efficiently is often more involved, especially with diffi- 
cult non-invertible densities with infinite support. Since rarely needed, variates 
from the tail can safely be generated by a slower method |U [57J 12H] ■ Overall, 
a significant speed-up can be achieved. In this paper we show how to sample 
directly and efficiently via a rejection technique a random number X such that 
X > x or X < x where x G (—00, +00) at least within the limits of the numerical 
representation. This is achieved by using properties of the transform representa- 
tion of the distributions. The examples we use for demonstration are the Levy 
a-stable [3TJ [3TJ [35] and the Mittag-LefHer one-parameter probability densi- 
ties [T5] . A transformation formula for the former is well known [2 [J2] , while 
the transform representation of the latter was discovered [5|. IT6 | IT7 | fl^ l rT9l . [20 | [33] 
and applied [SJ [HI HO] only recently. The two distributions are generalizations of 
the Gaussian and exponential distribution respectively and play an important 
role together for the solution of the space-time fractional diffusion equation. 

Our rejection concept is general for any distribution that provides a trans- 
form representation. It can sample efficiently from arbitrary infinite or finite 
intervals as opposed to other existing methods that are designed especially for 
certain densities. In this work we do not consider the technical details of a 
speed-optimized implementation, but explain the basis of the algorithm and 
show example applications. The method is based on properties of the two- 
dimensional transform maps that seem unnoticed yet. 

The assumption for using the method introduced here is that the tail re- 
gion requires high accuracy due to high demands on statistics as well as speed. 
The transformation formula by Chambers, Mallows and Stuck [5] for example 
is exact and for most applications the recommended method for the produc- 
tion of Levy a-stable random numbers [32]. The replacement of the tails by 
a simple invertible Pareto function is not totally appropriate because this is 
only an asymptotic approximation; moreover it introduces a transition region. 
The more sophisticated and smooth this transition, the more complicated and 
slower the overall procedure. Such a replacement of the tail contrasts the initial 
goal of speed. But the most demanding contemporary applications of random 
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numbers [551 S2], as of the two suitable examples we treat here, will require 
large amounts and therefore fast production. The tails should be accurate with- 
out an approximated transition region from the density body to its tails. In 
some cases fast series expansion methods can be used but with a compromise 
in accuracy [I]. A more detailed analysis of such considerations can be found 
in Ref. [39 , where most known algorithms for the Gaussian distribution (as a 
simple and special case of the Levy a-stable distribution) are analyzed in the 
context of contemporary statistical applications as well as expectations of future 
demands. It is argued extensively how speed of production implies the demand 
for very many random numbers, which in turn requires greater accuracy of the 
resulting distribution. 

Consider the Ziggurat rejection method by Marsaglia et al. US] [57] that 
was introduced to produce Gaussian and exponential random numbers. It is 
an exact method up to the numerical limits of floating point representation. 
In principle it is applicable to all decreasing or symmetric densities, provided a 
suitable tail sampling method is available [24 . In particular the implementation 
by Marsaglia and Tsang [2SJ and a recent version by Rubin and Johnson [35] are 
about two orders of magnitude faster on contemporary processors than other 
dedicated methods for Gaussian and exponential random variates; therefore it 
is likely to outrun any non-trivial transformation method by at least the same 
factor. The hurdles to apply the Ziggurat method to other densities with infinite 
support, with additional parameters and for which no closed form or simple 
transform exist are: a) the costly setup of the look-up table, b) the necessity 
of equal areas of the rectangles covering the density as well as the area under 
the tail and finally c) a reasonably fast and accurate tail sampling method. 
Difficulty a) must be evaluated in relation to the required number of variates if 
it is possible to predict the setup costs as a function of the density parameters. 
The meaning of "fast" in c) is defined by the ratio of tail variates versus body 
variates and the speed of the body sampling method. A slow tail sampling can 
always be balanced by sufficiently infrequent calls to the latter. 

Provided the complementary cumulative density function J f{x') dx' can 
be computed sufficiently exact on demand, then any required value of the tail 
surface, and thus any relative frequency of calls to the tail sampling function, can 
be achieved in the setup of the Ziggurat by an iterative process. For the details 
of the setup refer to Ref. [25] and for alternative concepts to Ref. [35] . Indepen- 
dently of such considerations the production of Levy a-stable random numbers 
in the tails, but also in arbitrary finite intervals, are themselves examples where 
the method introduced in this paper is suitable. Of course the Ziggurat method 
is applicable to non-symmetric decreasing densities by representing two halves 
with separate generators which have to be called alternatingly in a ratio that 
corresponds to the ratio of respective areas covered by each halves. 

In Sec. [2] we introduce the Levy a-stable probability density on the basis of 
which Sec. [3j explains our method. In Sec. [4] the Mittag-Leffler distribution, its 
transform representation and transform map are presented. 
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L af3j s(x) = — I 4> al3l s(k) exp(-ikx) dk, (1) 



2 The Levy a-stable probability density and its 
transform map 

A convenient representation of the Levy probability density function in its most 
popular parametrization [3ll 132J 02] is via the inverse Fourier transform of its 
characteristic function: 

.-fVi = 
where 

( - 7 Q |fc| Q (l - ^sign(fe) tan(f a)) + iSk for a ^ 1 
log (f> a 0fS(k) = < (2) 
[ — T| A; | (l + i/?sign(fc)f log \k\) + iSk for a = 1. 

The index or order a G (0,2] determines the exponent of the power-law tail. 
The parameter (3 G [—1,1] governs the skewness, 7 € (0,oo) the horizontal 
scale and 5 £ (—00,00) the location. The advantage of this parametrization is 
that the density and the distribution function are jointly continuous in all four 
parameters; the same applies to the convergence to the power-law tail. The last 
two parameters can safely be set to 1 and without loss of generality. Other 
values can be obtained through 

XafjjS = 7X Q/ 3io + S. (3) 

We therefore omit 7 and 6 in the subscripts and also (3 if equal to zero. The 
symmetric case with (3 — has the simpler form of an inverse cosine transfor- 
mation 

1 f°° 

L a (x) = - / exp(-fc Q ) cos(fcx) dfc. (4) 







Rejection methods for Levy a-stable random numbers that use asymptotic series 
representations of the density function are sometimes used if speed has highest 
priority [1] . However, the known types of series expansions for the Levy density 
tend to become inaccurate especially in the tails and also account for a certain 
fraction of uniform random numbers to be lost (rejected) in the sampling. To 
achieve best performance (minimum rejection rate and maximum accuracy) one 
must use different versions of the algorithms and expansions depending on the 
combination of parameter values and their range. This is the case in particular 
for (3 ^ 0. A review on these methods and their deficiencies can be found in 
Ref. g. 

A transformation method for Levy a-stable random numbers by Chambers, 
Mallows and Stuck has been available for over 30 years j2]. Two independent 
uniform random numbers {/, V G (0, 1) are mapped via a transform F a p(U, V) 
such that X = F a/3 (U,V) is distributed correctly according to L a p{x). The 
general case for a ^ 1 is given by 

X = F aP (U, V) = M°(* + *o)) ( -logt/cos* \ ^ 
at>y ' ; cos$ Vc° s (*-a(* + *o))/ 

where $ = 7r (V — |) and <&o = \ n P^~ L ^~^ L 1 while for a = 1 



4 



The symmetric case with (3 — simplifies to 



x = Fai u, v) = (- l T co :l) " 1/Q ■ W 

V ; cos$ VcosU 1 - / 
The variables X±, ...,Xn are stable as well as their normalized sum 

1 N 

i=i 

This transform representation is a mixture of the form g(V)W 1 ~ 1 ^ a where g(V) 
is a real valued random number and W is exponentially distributed. Figs. [2] 
and [3] show symmetric and asymmetric examples of the mapping of the random 
number plane (U, V) to "quantiles" of the probability density via the map X — 
F a p(U,V). Colors are used to designate the respective regions X\ < X < Xi + \ 
separated by isolines defined by dF a p(U, V) = 0. The pictures show isolincs 
as borders between colors for Xi — 0, ±0.5, ±1, ±1.5, ... . The colors in the 
map and in the respective histogram correspond to each other and all points 
(U, V) on the same isoline are mapped onto exactly one unique number. Fig. [4] 
shows the behaviour of the isolines only with further decreasing a. Notable is 
the analytic Cauchy case a = 1 whose inversion formula depends only on one 
variable. This is expressed by perfectly vertical isolines. For values of a < 1 
the overall behaviour turns over and the slopes change sign in each half of the 
unit square. The pictures showing isolines are produced with Matlab 7.4's 
contourf function on a grid of size 800x800. 

We would like to remark that different solutions are thinkable of how to sam- 
ple uniform random points in a specific region in the ([/, V^)-plane. A differential 
equation for the isolines can be obtained via the implicit function theorem by 
Ulisse Dini [B]: 

% dF(u,v) , dF(u,v) , 

= dF(u,v) = y —^du+ )r^dv ; (9) 

ou ov 



rearranging 



dv(u) = S dF(u,v) \ dF(u, v) 
du \ dv J du 



With an appropriate initial condition this differential equation defines the isoline 
v(u) in the coordinate square spanned by u, v. The alternative representation of 
u as a function of v is equally appropriate from the mathematical point of view, 
but is less convenient in this case for symmetry reasons. We skip additional 
considerations on singularities and limiting behaviour. For a — 2 and f3 = 
Eq. reduces to X = F 2 (U,V) = 2V~ log Usm(Tr(V - 1/2)), which is the 
Box-Muller method for Gaussian deviates with standard deviation a = y/2. 
The corresponding map is shown in the upper left of Fig. [2| T he value x in the 
condition X > x determines the initial condition for Eq. (|10| that determines 
the isoline, i.e. x in the condition X > x, and for a = 2 it can be chosen 
on the boundary of the square U,V £ (0,1). Two other analytic limit cases 
for (3 = 0, where L a (x) can be written in terms of elementary functions, are 
the Cauchy distribution, with a = 1 and X = F\(U) = tan(7r(J7 — 1/2)), 
and the Levy distribution, with a = 1/2 and X = F 1 / 2 (U,V) = — tan(7r(V r — 
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Figure 2: (Color online) The map X = F a (U,V) with U,V € (0,1) giving 
the symmetric Levy distribution L a (x) for different values of a. For a = 2 
the picture corresponds to the Box-Muller map for the generation of Gaussian 
random numbers. The bottom part of each map shows the respective histogram. 
Areas with equal colors correspond to each other. Note that the transition from 
a = 2 to a < 2 is discontinuous for u = and u = 1 and the points (0,1) and 
(1,1) develop a singularity. 
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Figure 3: (Color online) The map X - 
distribution L a p(x) for two values of [3 
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F a p(v,,v) giving the asymmetric Levy 



l/2))/(21ogJ7cos(7r(^ - 1/2))). Note that for values of a ^ 2 the map F 
is singular in the points (0,1) and (1,1). In such cases the initial condition 
cannot be chosen on the boundary, which considerably complicates the situation 
numerically 



Starting from the simplest case, insertion of ^(u, v) into Eq. (10 1 yields 



dv{u) cot(irv(u)) 
du 2iru log(w) 



Insertion of F a (u, v) into Eq. (10 1 yields 



(11) 



dv(u) 
du 



(a -I) 



TTU log(u) 



tan(^( V ( U )- -)(!-«)) (a -l) 2 



+ cot(7n;(u)) — a cot 7r v(u) — 



(12) 



One way to sample directly and uniformly from the area under v(u) would be 
an area-preserving map of a square domain spanned by two uniform random 
numbers, e.g. U, V £ (0, 1), or any other suitable two dimensional domain onto 
this area. To our knowledge this solution is not available yet. Alternatively, 
the function v(u) can be obtained numerically via integration or by appropriate 
algorithms for the generation of isolines. Once data points for v(u) are obtained, 
any method that samples uniformly the region X < x or X > x is suitable 
in principle. With this, the generation of a tail variable constitutes in itself 
a standard non-uniform variate generation task. It is the initial scenario of 
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Figure 4: Isolines of the map X = F a (u,v) with u,v e (0,1) for decreasing 
values of a. The regions with increasingly divergent gradient (upper corners) is 
not shown beyond |a;| > 600. Note that the orientation of the isolines flips over 
with decreasing values of a at exactly a = 1. 
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Figure 5: (Color online) Intuitive, coarsely tiled, example of the tiling in the 
u-v square for sampling symmetric Levy a-stable random variates with the 
condition X = F(u, v) < — 1 and a — 1.8, (3 — 0. The tiled area can be sampled 
efficiently while only points in the red shaded region are rejected. Tiles with 
direct acceptance do not require the acceptance comparison X — F(u, v) < — 1. 



sampling uniformly under a curve, but with the great simplification of a finite 
support. However, this is not the route we propose for three reasons. First, the 
numerical solution of Eq. (12 1 is cumbersome. Second, the initial condition has 
to be found within the u-v square due to the above mentioned singularities. The 
subsequent integration in two directions must be guaranteed to work unattended 
and automatically as a black box with a and (3 as the only parameters. Third, 
the outcome is not exact in the sense that the sampled random tail variates are 
distributed with respect to an approximated probability density function based 
on the data point representation of the isoline. As it will turn out a numerical or 
analytic representation of the isoline is not a required piece of information and its 
calculation can be avoided. It can also be shown that the isolines are monotonic 
in u in the regions F a p(u, v) < and F a p{u, v) > which is a useful property in 
Sec. [3] Although the approximation of density functions is commonly accepted 
as a reasonable compromise in several applications we introduce in the next 
section a simple graphical method without this disadvantage. 



3 Sampling method and example application 

We introduce the method using simple intuitive examples. The production 
algorithm relies on the rejection method whose invention dates back to von 
Neumann |40] and which we do not rehearse here. Fig. [5] demonstrates a com- 
putationally efficient concept for uniform sampling in a certain two-dimensional 
region. In the first example we aim at producing Levy a-stable random vari- 
ates with parameters a = 1.8, (3 = and the condition X < —1. The map 
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Figure 6: (Color online) Two different tiling refinements of the region corre- 
sponding to the condition X < — 12 which is a narrow strip along the left and 
bottom of the unit square. Only the lower left corner is shown on a scale that 
magnifies the tiling to a visible size. The row of tiles on the bottom samples 
a narrow strip below the isoline. In the right panel the rejection rate is sig- 
nificantly lower. It may help intuition that the colored dots are the uniformly 
distributed pairs (u, v). 

Fi.$(u,v) for this choice of parameters is also shown in Fig. [2j It corresponds 
to a relatively large region in the left part of the square. 

We performed a straightforward and simple tiling of this region using square 
tiles that can be refined, for example, iteratively maintaining complete cover- 
age while minimizing the excess area of the tiles that stick out of the region 
defined by F±,s(u, v) < x. Uniform sampling of the tiled area accepting all 
X = Fi.s(u,v) < x and rejection of all other samples achieves the desired tail 
sampling. The size of the tiles can be chosen to achieve an arbitrarily low rejec- 
tion rate. In the example shown in Fig. [5] the tiling is refined only moderately 
to convey the situation. For tiles that lie completely underneath the isoline the 
test Fi,$(u,v) < x must not be executed. With dense tiling this comparison is 
therefore hardly needed and indeed must be avoided to yield a speed-up with re- 
spect to the transformation method. The setup and production loop of random 
variates follows the steps: 



10 



Algorithm Input: 16R. 
0: Setup: 

Tiling of the region F(u, v) < x using a method of choice. 
Label tiles with an integer index. 

Label tiles that are intersected by the isoline F(u, v) = x. 

1: Draw a random integer tile index with uniform probability. 

2: Draw a random coordinate (u, v) with uniform probability within this tile. 

3: Test if the tile is intersected by the isoline (table look-up). 

If yes, go to 4. If no, accept X — F(u, v) and go to 1 (direct acceptance). 

4: Test if (it, v) satisfies X — F(u, v) < x. 

If yes, accept X, otherwise reject and go to 1. 
Note that for monotonic isolines the position of a tile with respect to an isoline, 
i.e. whether underneath, above or intersected, can be determined by evaluating 
the map for at most two corners. Step 4 is unlikely to be carried out if the 
coverage is dense, giving nearly a zero rejection rate. Overall, this procedure 
is efficient in setup and production for a sufficiently dense tiling. Furthermore, 
with small modifications of the above acceptance and rejection conditions in the 
pseudo code, the tiling and production of random numbers on a finite interval 
X G [0:1,0:2] is geometrically and algorithmically equivalent to generating num- 
bers from the tail. This requires the tiling of a region in the u-v square between 
two isolines with the condition x\ < X < X2- 

Fig. [6] shows the map for the left tail regions of the u-v square with X < 
— 12, which is more realistic for the purpose of tail sampling. This condition 
corresponds to sampling a narrow strip at the bottom and left sides of the unit 
square. The figure only shows the corner at the origin. The bottom strip of 
tiles samples an extremely narrow strip than is not visible on this scale. The 
iterative tile refinement in the setup stage is acceptably fast, below a second in 
our non-optimized code, down to the level on the right panel of Fig.[6]to achieve 
a rejection rate below 1%. Different values of a > 0.1 as well as not too extreme 
values of x have no significant influence on the setup performance achieving 
a rejection rate of 1%. Note that the speed of random number production 
is independent of the number of tiles. In our case it amounts to 2.3 million 
tail variates per second on a PC with a 2.4 GHz Intel Pentium 4 processor 
using the GNU C++ compiler version 3.2.2 Linux and optimization level -03. 
As the uniform random number generator we used the XOR shift SHR3 by 
Marsaglia [55]. The colouring of the acceptance and rejection regions in Figs. [5] 
and [6] are produced by green and red coloured dots representing the random 
uniform coordinates (u, v). 

We would like to stress that the method of tiling as well as the form of the 
tiles is in principle arbitrary. Equal size and shape is computationally advanta- 
geous, but this issue is not the focus of the present work. Of course any tiling 
technique that produces a similar result is suitable, using either square or rect- 
angular tiles. However, the choice of square equal tiles is algorithmically very 
simple and likely to outrun an adaptive scheme with more complex shapes in 
setup and also production. The iterative tiling refinement, as performed in the 
above examples, is robust and fast also for large values of \x\. The rejection 
scheme is in principle similar to the Ziggurat implementation of Ref. [28 . It 
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also needs the setup of a data structure that covers a region by equal area rect- 
angles. The details of the tiling method that is more general and applicable to 
random number production directly via the probability density are described in 
separate work [7]. 



4 The Mittag-Leffler probability distribution 

Our second example density is less know in scientific applications, even less so 
its transform. The Mittag-Leffler probability distribution appears e.g. in the 
analytic solution of the time-fractional Fokker-Planck equation j8j [T2l [13j [14] . 
The generalized Mittag-Leffler function is defined as jTTJ [T5] 

oo n 

For our purposes it is sufficient to restrict the example to the one-parameter 
Mittag-Leffler function which plays an important role in the stochastic solution 
of the time-fractional diffusion equation. The series representation is 

E °W = £f(£+iy zeC ' (14) 

n— y 

leading to a pointwise representation on a finite interval. The Mittag-Leffler 
function with argument z — —t a , t € K reduces to a standard exponential 
decay, e - *, with a = 1; when < a < 1, the Mittag-Leffler function is approxi- 
mated for small values of t by a stretched exponential decay (Weibull function) , 
exp(— t a /T(l + a)) and for large values of t by a power law, t~ a /T(l — a), 
see Fig. [jj top left plot. There is increasing evidence for physical phenomena 
[3 [501 EH ELI] and human activities [TJ [35J [37] that do not follow neither expo- 
nential nor, equivalently, Poissonian statistics. The Mittag-Leffler distribution 
is an example of power-law distributed waiting times. They arise as the natural 
survival probability leading to time-fractional diffusion equations. 



Eq. (14 1 is the complementary cumulative distribution function, also called 



survival function; the proability density is 

- ^E a {-n. (is) 

In past applications Mittag-Leffler random numbers were produced by rejection 



through a pointwise representation via Eq. (14 1, which is inefficient due to the 
slow convergence of the series. In some cases concepts to avoid Mittag-Leffler 
random numbers were presented 22, 23 due to the difficulty of their production. 
In this context it had not been recognized immediately that transformation 
formulas analogous to Eq. ([jj are available [H [HI CH LIZ1 HH1 HSJ HHJ I3H]- The 
most convenient expression is due to Kozubowski and Rachev |20j : 

T = M a (U, V) = log U f f Sm , (mr ) - cos(aTr)) ' , (16) 
\tan(airV ) J 

where U, V £ (0, 1) are independent uniform random numbers and T is a Mittag- 



Leffler random number. For a = 1, Eq. (16 1 reduces to the transform for the 
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exponential distribution: Mi (U,V) = log U. Fig. [7] shows the map M a (U, V) of 
the transform representation Eq. (16 1 as borders between intervals correspond- 
ing to t = 0, ±0.5, ±1, ±1.5, ... The exponential case with a = 1 depends on only 
one random variable in the u-v square which is expressed by perfectly horizontal 
isolines. For a < 1 the left and right edges develop singularities. It is not recom- 
mended to use Eq. (14) and summation of many terms for the computation of 
E a (—t a ). A more elegant and accurate method is presented in Ref. [TTl IMl ITS] . 
For the generation of random numbers we use the implementation in Ref. [10 . 
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Figure 7: The Mittag-LefHer function E a (—t a ) in a log-log plot (top left) and 
the transformation map X = M a (U,V), Eq. (|T6]) , in terms of isolines for five 
values of a. The case a = 1 corresponds to the standard exponential function. 
The regions on the left side of the maps is not shown beyond t > 600 due to an 
increasingly divergent gradient. The two plots at the bottom repeat the case 
with a = 0.9 using colors and showing the corresponding histogram. 
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5 Summary and conclusion 



We have demonstrated some properties of the Chambers-Mallows-Stuck and 
Kozubowski-Rachcv transformation maps exemplifying the production of ran- 
dom numbers with the former. The interpretation as a two-dimensional map 
from the unit square to the real numbers allows to associate arbitrary intervals 
on the support of the density with well defined finite regions of the map do- 
main. The uniform sampling of such regions produces directly random numbers 
exactly within the respective intervals. We have also introduced an efficient 
concept for the automatic setup of a random number generator that makes use 
of this property. The resulting generator can in principle produce random num- 
bers in intervals which can be disconnected and any combination of the kind 
(—00, xi] U [X2, X3] U . . . U [x n , +00), Xi € R. Most importantly, the sampling of 
tails as shown here can be used as a tail handling method in fast implementations 
of random number generators for which a candidate is the Ziggurat method im- 
plementation that was proven to greatly outrun simple inversion methods. The 
present work opens the route for the speedup of many known random number 
generators that rely on transform representations. 
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